{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real Test of the nonreciprocal RCWA\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(220,)\n"
     ]
    }
   ],
   "source": [
    "## MATERIALS\n",
    "def Lorentz_SiC(wvlen):\n",
    "\n",
    "    # lambda in micron\n",
    "    mu_0=4*pi*1e-7;  #N/A^2\n",
    "    e_0=8.854e-12;\n",
    "\n",
    "    wn = 1e4/wvlen; # wavenumber [cm-1]\n",
    "    omega_L=969; #[cm-1]\n",
    "    omega_T=793; #[cm-1]\n",
    "    Gamma=4.76;  #[cm-1]\n",
    "    ep_inf=6.7;\n",
    "    ep = ep_inf*(1+(omega_L**2-omega_T**2)/(omega_T**2-wn**2-1j*Gamma*wn));\n",
    "\n",
    "    n = np.real(np.sqrt(ep));\n",
    "    k = np.imag(np.sqrt(ep));\n",
    "    return n,k;\n",
    "\n",
    "def MOInAs(wvlen,B):\n",
    "    #eps_yy = eps_xx\n",
    "    #lambda in microns\n",
    "    #B in tesla\n",
    "\n",
    "    lightv=2.99792458e8;\n",
    "    hbar=1.054571726e-34;\n",
    "    kb=1.3806488e-23;\n",
    "    mu0=4*np.pi*1e-7;\n",
    "    eps0=8.854187817e-12;\n",
    "    me =9.10938291e-31;\n",
    "    eV =1.60217657e-19;\n",
    "\n",
    "    eps_inf=12.37;\n",
    "\n",
    "    c0=2.99792458e8;\n",
    "    w=2*np.pi*c0/(wvlen*1e-6);\n",
    "    n=7.8e17*1e6;\n",
    "    m=0.033*me;\n",
    "    alpha_factor=0.038197513;\n",
    "\n",
    "    wmega_p=np.sqrt(eV**2*n/m/eps0);\n",
    "    wmega_c=eV*B/m;\n",
    "\n",
    "    #single wavelength situation\n",
    "    lambda_list=wvlen*1e-6;\n",
    "\n",
    "    # perform the constants calculation\n",
    "    wmega_list = 2*np.pi*lightv/lambda_list;\n",
    "    alpha_exp=alpha_factor*100*(lambda_list*1e6)**3;\n",
    "    # this determination is very accurate already, even near the ENZ region\n",
    "    tau= wmega_p**2/(wmega_list**2*lightv*alpha_exp*np.sqrt(eps_inf-wmega_p**2/wmega_list**2 +0*1j));\n",
    "\n",
    "    denom= wmega_list*( (1j+wmega_list**tau)**2-(wmega_c*tau)**2 );\n",
    "\n",
    "    eps_xx= (eps_inf-wmega_p**2*tau*(1j+wmega_list*tau)/denom          );\n",
    "    eps_xy= (     1j*wmega_p**2*tau*(wmega_c*tau)      /denom          );\n",
    "    eps_zz= (eps_inf-wmega_p**2.*tau/( wmega_list*(1j+wmega_list*tau) ) );\n",
    "\n",
    "    return eps_xx, eps_xy, eps_zz\n",
    "\n",
    "wvlen_vec = np.concatenate((np.linspace(24.2,24.59,20),np.linspace(24.6, 25.1, 180), np.linspace(25.11, 25.5, 20)));\n",
    "N_wvlen = len(wvlen_vec);\n",
    "print(wvlen_vec.shape)\n",
    "\n",
    "Bfield = 0.1;\n",
    "\n",
    "for ind in range(N_wvlen):\n",
    "    #print(wvlen_vec[ind])\n",
    "    [eps_mo_11, eps_mo_12, eps_mo_33] = MOInAs(wvlen_vec[ind],Bfield);\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structure specification\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-17-acb27d565db8>, line 4)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-17-acb27d565db8>\"\u001b[1;36m, line \u001b[1;32m4\u001b[0m\n\u001b[1;33m    eps0 = 8.854e-12 * L0;  % vacuum permittivity in farad/L0\u001b[0m\n\u001b[1;37m                            ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "theta = 65 * np.pi/180;\n",
    "\n",
    "##\n",
    "eps0 = 8.854e-12 * L0;  % vacuum permittivity in farad/L0\n",
    "mu0 = pi * 4e-7 * L0;  % vacuum permeability in henry/L0\n",
    "c0 = 1/sqrt(eps0*mu0);  % speed of light in vacuum in L0/sec\n",
    "# omega = 2*pi*c0/wvlen;  % angular frequency in rad/sec\n",
    "\n",
    "Bfield=0.1;\n",
    "\n",
    "## Specify the material geometry (can be done outside for loop)\n",
    "d_sub = 1.001;\n",
    "d_mo = 20;\n",
    "d_SiC = 4;\n",
    "d_grating = 0.247;\n",
    "w_grating = periodl*0.47+1e-8;\n",
    "\n",
    "within_sub = @(x, y) y <= d_sub;\n",
    "within_mo = @(x, y) y > d_sub & y <= d_sub + d_mo;\n",
    "within_SiC = @(x, y) y > d_sub + d_mo & y <= d_sub + d_mo + d_SiC;\n",
    "within_grating = @(x, y) y > d_sub + d_mo + d_SiC & y < d_sub + d_mo + d_SiC + d_grating & abs(x) < w_grating/2;\n",
    "\n",
    "# %% Solve for the reflection at each frequency\n",
    "# % N_wvlen = 32*4;\n",
    "# % wvlen_vec = linspace(24.86, 25.02, N_wvlen);\n",
    "\n",
    "wvlen_vec = [linspace(24.2,24.59,20),linspace(24.6, 25.1, 180), linspace(25.11, 25.5, 20)];\n",
    "N_wvlen = length(wvlen_vec);\n",
    "R_vec = [];\n",
    "Bfield = 1;\n",
    "\n",
    "for ind = 1:N_wvlen\n",
    "    disp(ind);\n",
    "    [eps_mo_11, eps_mo_12, eps_mo_33] = MOInAs(wvlen_vec(ind),Bfield);\n",
    "\n",
    "    ## substrate material PEC, or 'silver'\n",
    "    eps_Ag = -1i*1e8\n",
    "    \n",
    "    ## silicon carbide (grating material)\n",
    "    eps_SiC(ind) = np.conj(Lorentz_SiC(wvlen_vec[ind]));\n",
    "\n",
    "    \n",
    "    #run anisotropic simulation\n",
    "    \n",
    "    #store the refelection\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
